/** Runs Monte Carlo comparison of KK and SUR.
 ** Jason Wittenberg, 12/20/2001 witty@polisci.wisc.edu
 **/

/*** Set up Gauss ***/
new;
library cml,lpack,comp2,probs;

#include gauss.ext; 
#include cml.ext;
#include comp2.ext;

cmlset;

rndseed 123456789;            /* I added this-- fix seed */
setprior=0;                   /* 1=Use KK's prior, otherwise 0 */

/** Set Estimation Options **/

_IncRun_Verbose = 0;

v1hatKK={}; 
v2hatKK={}; 
v1hatSUR={}; 
v2hatSUR={};

b10kk={};
b11kk={};
b12kk={};
b20kk={};
b21kk={};
b22kk={};

b10sur={};
b11sur={};
b12sur={};
b20sur={};
b21sur={};
b22sur={};

ndist = 500;                             @ number of districts @ 
nmonte = 1000;                           @ number of monte carlos @ 
p=2;

x1 = rndmn(1,1,ndist);                  @ explanatory var 1 @ 
x2 = rndmn(1,1,ndist);                  @ explanatory var 2 @ 
X = ones(ndist,1)~x1~x2;                @ matrix of X's, including constant @ 
mu1 = .2 + .2*x1 + .5*x2;               @ E(y1), the mean for party 1 @ 
mu2 = .8 - .7*x1 + .5*x2;               @ E(y2), the mean for party 2 @ 
v1 = exp(mu1)./(exp(mu1)+exp(mu2)+1);   @ "true" vote share for party 1 @ 
v2 = exp(mu2)./(exp(mu1)+exp(mu2)+1);   @ "true" vote share for party 2 @ 
v3 = 1 - v1 - v2;                       @ "true" vote share for party 3 @ 
nu = 5;                                 @ degrees of freedom, as in HKK @ 
sigma = (.06~.06)|(.06~.08);

"Initial Time: " time;
  
for i (1,nmonte,1); 
                                               @ Step 1: Prepare a dataset @ 
   y1={}; 
   y2={};    
   for j (1,ndist,1); 
      e = rndmt(zeros(2,1),sigma,nu,1);        @ epsilons for the 2 equations @ 
      y1 = y1 | (mu1[j] + e[1]);         
      y2 = y2 | (mu2[j] + e[2]); 
   endfor; 
   vv1 = exp(y1)./(exp(y1)+exp(y2)+1);
   vv2 = exp(y2)./(exp(y1)+exp(y2)+1);
   
                                               @ Step 2: Do KK @
   dta=0;
   
   dta = lput(dta,vv1~X,1);
   dta = lput(dta,vv2~X,2);
   
   _cml_ParNames = { Const_1, X1_1, X2_1, Const_2, X1_2, X2_2, Tdf, TS11, TS12, TS22 };
    
   { theta,mlik,grad,cov,retcode } = RegMT2(dta);                @ From KK replication DS @
   
   Const_1=theta[1];                                             @ Collect KK effect params @
   X1_1=theta[2];
   X2_1=theta[3];
   Const_2=theta[4];
   X1_2=theta[5];
   X2_2=theta[6];
   
   Tdf = tdfi(p,theta);                                          @ Degrees of freedom @
   Sig = tsigi(p,theta);
   
   ts11=sig[1,1];                                                @ Elements of sigma @
   ts22=sig[2,2];
   ts12=sig[1,2];
   
   yhat1 = Const_1 + X1_1*x1 + X2_1*x2;                          @ dimension ndist x 1 @ 
   yhat2 = Const_2 + X1_2*x1 + X2_2*x2;                          @ dimension ndist x 1 @ 
   
   b10kk=b10kk | Const_1;                                        @ Betas from each run @
   b11kk=b11kk | X1_1;
   b12kk=b12kk | X2_1;

   b20kk=b20kk | Const_2;
   b21kk=b21kk | X1_2;
   b22kk=b22kk | X2_2;
   
   v1hatKK = v1hatKK ~ (exp(yhat1)./(exp(yhat1)+exp(yhat2)+1));  @ ndist x ... @ 
   v2hatKK = v2hatKK ~ (exp(yhat2)./(exp(yhat1)+exp(yhat2)+1));  @ ndist x ... @ 
   
                                                @ Step 3: Do SUR @
   _output = 0;                                                  @ suppress OLS output @
                                                    
   {vnam,m,b1,stb,vc,stderr,sigma,cx,rsq,resid,dwstat} = ols(0,y1,X); @ OLS on party 1 @
   {vnam,m,b2,stb,vc,stderr,sigma,cx,rsq,resid,dwstat} = ols(0,y2,X); @ OLS on party 2 @
   
   yhat1 = b1[1] + b1[2]*x1 + b1[3]*x2;                          @ dimension ndist x 1 @ 
   yhat2 = b2[1] + b2[2]*x1 + b2[3]*x2;                          @ dimension ndist x 1 @ 
   
   
   v1hatSUR = v1hatSUR ~ (exp(yhat1)./(exp(yhat1)+exp(yhat2)+1));  @ ndist x ... @ 
   v2hatSUR = v2hatSUR ~ (exp(yhat2)./(exp(yhat1)+exp(yhat2)+1));  @ ndist x ... @ 
   
   b10sur=b10sur | b1[1];
   b11sur=b11sur | b1[2];
   b12sur=b12sur | b1[3];

   b20sur=b20sur | b2[1];
   b21sur=b21sur | b2[2];
   b22sur=b22sur | b2[3];

    
endfor; 

v3hatKK = ones(ndist,nmonte) - v1hatKK - v2hatKK; 
v3hatSUR = ones(ndist,nmonte) - v1hatSUR - v2hatSUR; 

v1=100.*v1;                                                      @ Transform votes to % @
v2=100.*v2;
v3=100.*v3;
v1hatkk=100.*v1hatkk;
v2hatkk=100.*v2hatkk;
v3hatkk=100.*v3hatkk;
v1hatsur=100.*v1hatsur;
v2hatsur=100.*v2hatsur;
v3hatsur=100.*v3hatsur;


v1biasKK = meanc((v1hatKK - v1)');                               @ kk bias for each district @
v2biasKK = meanc((v2hatKK - v2)');  
v3biasKK = meanc((v3hatKK - v3)');
v1varKK = stdc(v1hatKK')^2;                                      @ kk variance for each district @
v2varKK = stdc(v2hatKK')^2;
v3varKK = stdc(v3hatKK')^2;
v1mseKK = v1biasKK^2 + v1varKK;                                  @ kk mse for each district @
v2mseKK = v2biasKK^2 + v2varKK;
v3mseKK = v3biasKK^2 + v3varKK;

v1biasSUR = meanc((v1hatSUR - v1)');                             @ sur bias for each district @
v2biasSUR = meanc((v2hatSUR - v2)'); 
v3biasSUR = meanc((v3hatSUR - v3)'); 
v1varSUR = stdc(v1hatSUR')^2;                                    @ sur variance for each district @
v2varSUR = stdc(v2hatSUR')^2;
v3varSUR = stdc(v3hatSUR')^2;
v1mseSUR = v1biasSUR^2 + v1varSUR;                               @ sur mse for each district @
v2mseSUR = v2biasSUR^2 + v2varSUR;
v3mseSUR = v3biasSUR^2 + v3varSUR;

/* Now save effect parameters--- are they biased? */

b10biaskk=meanc(b10kk-0.2);
b10biassur=meanc(b10sur-0.2);

b11biaskk=meanc(b11kk-0.2);
b11biassur=meanc(b11sur-0.2);

b12biaskk=meanc(b12kk-0.5);
b12biassur=meanc(b12sur-0.5);

b20biaskk=meanc(b20kk-0.8);
b20biassur=meanc(b20sur-0.8);

b21biaskk=meanc(b21kk+0.7);
b21biassur=meanc(b21sur+0.7);

b22biaskk=meanc(b22kk-0.5);
b22biassur=meanc(b22sur-0.5);

b10varkk=stdc(b10kk)^2;
b11varkk=stdc(b11kk)^2;
b12varkk=stdc(b12kk)^2;
b20varkk=stdc(b20kk)^2;
b21varkk=stdc(b21kk)^2;
b22varkk=stdc(b22kk)^2;

b10varsur=stdc(b10sur)^2;
b11varsur=stdc(b11sur)^2;
b12varsur=stdc(b12sur)^2;
b20varsur=stdc(b20sur)^2;
b21varsur=stdc(b21sur)^2;
b22varsur=stdc(b22sur)^2;

b10msekk=b10biaskk^2+b10varkk;
b11msekk=b11biaskk^2+b11varkk;
b12msekk=b12biaskk^2+b12varkk;
b20msekk=b20biaskk^2+b20varkk;
b21msekk=b21biaskk^2+b21varkk;
b22msekk=b22biaskk^2+b22varkk;

b10msesur=b10biassur^2+b10varsur;
b11msesur=b11biassur^2+b11varsur;
b12msesur=b12biassur^2+b12varsur;
b20msesur=b20biassur^2+b20varsur;
b21msesur=b21biassur^2+b21varsur;
b22msesur=b22biassur^2+b22varsur;

"Bias for party 1";
" KK    SUR ";
" ";
b10biaskk~b10biassur;
b11biaskk~b11biassur;
b12biaskk~b12biassur;

"Bias for party 2";
" KK    SUR ";
" ";
b20biaskk~b20biassur;
b21biaskk~b21biassur;
b22biaskk~b22biassur;


"Final time: " time;

fullmc=v1~v2~v3~v1hatkk~v2hatkk~v3hatkk~v1hatsur~v2hatsur~v3hatsur;
save fullmc;                          @ Save to reproduce Figure 1 @

ds=v1biaskk~v2biaskk~v3biaskk~v1varkk~v2varkk~v3varkk~v1msekk~v2msekk~v3msekk~
v1biassur~v2biassur~v3biassur~v1varsur~v2varsur~v3varsur~v1msesur~v2msesur~v3msesur;


/* The following section computes the auxiliary (GK) quantities of interest */

surpld1 = {};                         @ Probability of losing deposit @
surpld2 = {}; 
surpld3 = {}; 
kkpld1 = {}; 
kkpld2 = {}; 
kkpld3 = {}; 

pcpsur={};                            @ Percent correctly predicted @
pcpkk={};
for ce (1,nmonte,1);
    tempsur=maxindc((v1hatsur[.,ce]~v2hatsur[.,ce]~v3hatsur[.,ce])');
    pcpsur=pcpsur~tempsur;
    tempkk=maxindc((v1hatkk[.,ce]~v2hatkk[.,ce]~v3hatkk[.,ce])');
    pcpkk=pcpkk~tempkk;
endfor;

  pcptempsur=pcpsur;                  @ N x MSIMS--each the winner of the simulated election @
  pcptempkk=pcpkk;
  
  truth=maxindc((v1~v2~v3)');         @ N x 1 @
  
  pcpsur = (truth.==pcptempsur);      @ NxMSIMS vector of 1's and 0's vector. 1's=correctly predicted @
  pcpkk =  (truth.==pcptempkk); 
  
  pcpsur = sumc(pcpsur')./nmonte;     @ Sum up 1's and divide by montes to get frac for each dist @
  pcpkk = sumc(pcpkk')./nmonte; 

v1hatsur=v1hatsur';
v2hatsur=v2hatsur';
v3hatsur=v3hatsur';
v1hatkk=v1hatkk';
v2hatkk=v2hatkk';
v3hatkk=v3hatkk';

/* Compute probability of losing electoral deposit--FC  districts */ 

for ce (1,ndist,1); 
   surpld1=surpld1 | counts(v1hatsur[.,ce],4.99)./nmonte; 
   surpld2=surpld2 | counts(v2hatsur[.,ce],4.99)./nmonte; 
   surpld3=surpld3 | counts(v3hatsur[.,ce],4.99)./nmonte; 
    kkpld1= kkpld1 | counts(v1hatkk[.,ce],4.99)./nmonte; 
    kkpld2= kkpld2 | counts(v2hatkk[.,ce],4.99)./nmonte; 
    kkpld3= kkpld3 | counts(v3hatkk[.,ce],4.99)./nmonte; 
endfor; 
surpld=surpld1~surpld2~surpld3;
kkpld=kkpld1~kkpld2~kkpld3;

/* Compute probability of a candidate winning the election--FC  districts */ 

temp1sur={};
temp1kk={}; 
for ce (1,ndist,1); 
  tempsur=maxindc((v1hatsur[.,ce]~v2hatsur[.,ce]~v3hatsur[.,ce])');    @ MSIMSx1 @ 
  temp1sur=temp1sur~counts(tempsur,1|2|3)./rows(tempsur);              @ 3xN @ 
  tempkk=maxindc((v1hatkk[.,ce]~v2hatkk[.,ce]~v3hatkk[.,ce])');        @ MSIMSx1 @ 
  temp1kk=temp1kk~counts(tempkk,1|2|3)./rows(tempkk);                  @ 3xN @ 
endfor; 
pcwdkk=temp1kk';                          @ Prob candidate winning district--Nx3 @
pcwdsur=temp1sur';

/* Compute prob that district is competitive: all betw 30 and 40%--FC only */ 

competitivesur={};
competitivekk={}; 
for ce (1,ndist,1); 
  tempsur=(v1hatsur[.,ce]~v2hatsur[.,ce]~v3hatsur[.,ce])'; 
  tempkk=(v1hatkk[.,ce]~v2hatkk[.,ce]~v3hatkk[.,ce])'; 
  temp1sur=((tempsur[1,.].>=30 .and tempsur[1,.].<40) .and 
                    (tempsur[2,.].>=30 .and tempsur[2,.].<40) .and 
                    (tempsur[3,.].>=30 .and tempsur[3,.].<40));        @  1xMSIMS @
  temp1sur=sumc(temp1sur')/cols(temp1sur);                             @  1x1 for district i @
  competitivesur=competitivesur | temp1sur;  
  temp1kk=((tempkk[1,.].>=30 .and  tempkk[1,.].<40) .and 
                    (tempkk[2,.].>=30 .and tempkk[2,.].<40) .and 
                    (tempkk[3,.].>=30 .and tempkk[3,.].<40));          @ 1xMSIMS @ 
  temp1kk=sumc(temp1kk')/cols(temp1kk);                                @ 1x1 for district i @ 
  competitivekk=competitivekk | temp1kk; 
endfor;

output file=monte.out reset;
ds~v1~v2~v3~surpld~kkpld~pcwdkk~pcwdsur~competitivekk~competitivesur~pcpsur~pcpkk;
output off;
